{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Poisson-type Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`x = mg(A,b,elem)` attempts to solve the system of linear equations `Ax = b` for `x` using geometric multigrid solvers. Inside mg, a coarsening algorithm is applied. The method is designed for the system from several finite element descritzations of elliptic equations on a grid whose topology is given by the array `elem`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PDE\n",
    "The classic formulation of the Poisson equation reads as\n",
    "\n",
    "$$ - \\Delta u = f  \\text{ in }  \\Omega, \\qquad u  = g_D  \\text{ on }\n",
    "\\Gamma _D,  \\qquad  \\nabla u\\cdot n = g_N  \\text{ on } \\Gamma _N, $$\n",
    "\n",
    "where $\\partial \\Omega = \\Gamma _D\\cup \\Gamma _N$ and $\\Gamma _D\\cap \\Gamma _N=\\emptyset$. \n",
    "We assume $\\Gamma _D$ is closed and $\\Gamma _N$ open. The corresponding bilinear from $$ \n",
    "a(u,v) := \\int _{\\Omega} \\nabla u\\cdot \\nabla v\\, {\\rm dxdy}$$ will lead to a symmetric and positive definite (SPD) matrix. \n",
    "\n",
    "Variable diffusion coefficient is allowed, i.e. the operator $-\\nabla \\cdot( d \\nabla u)$. When $d$ is highly oscillatory or anisotropic, the resulting SPD matrix is hard to solve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Elements\n",
    "\n",
    "`x = mg(A,b,elem)` works for linear finite element and mesh `elem` is obtained by `uniformrefine`, `bisect`, or `uniformrefine3` by default. For other finite elements, more mesh structure, e.g., `edge` can be provided as `varargin`. If no extra input information, `mg` will try to guess the type of element by comparing the size of the system with the number of nodes and elements and construct needed data structure.\n",
    "\n",
    "For 3-D adaptive grids obtained by `bisect3`, `HB` is needed for the coarsening, and should be listed as the first parameter in varargin. Interesting enugh, for 2-D adaptive grids obtained by `bisect`, only `elem` is enough for the auotmatical coarsening. See [coarsen](../afem/coarsendoc.html).\n",
    "\n",
    "Here is a list of possible elements.\n",
    "\n",
    "    - mg(A,b,elem)                      work in most scenario \n",
    "    - mg(A,b,elem,option,edge)          2-D quadratic P2 element or linear CR element\n",
    "    - mg(A,b,elem,option,HB)            3-D linear P1 element on adaptive meshes\n",
    "    - mg(A,b,elem,option,HB,edge)       3-D quadratic P2 element on adaptive meshes\n",
    "    - mg(A,b,elem,option,HB,face)       3-D non-conforming CR element on adaptive meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More input and output\n",
    "\n",
    "`x = mg(A,b,elem,options)` specifies options in the following list.\n",
    "\n",
    "    - option.x0: the initial guess. Default setting x0 = 0.\n",
    "    - option.tol: the tolerance of the convergence. Default setting 1e-8.\n",
    "    - option.maxIt: the maximum number of iterations. Default setting 200.\n",
    "    - option.N0: the size of the coarest grid. Default setting 500.\n",
    "    - option.mu: smoothing steps. Default setting 1.\n",
    "    - option.coarsegridsolver: solver used in the coarest grid. Default\n",
    "      setting: direct solver.\n",
    "    - option.freeDof: free d.o.f\n",
    "    - option.solver: various cycles and Krylov space methods\n",
    "        * 'NO'     only setup the transfer matrix\n",
    "        * 'Vcycle'      V-cycle MultiGrid Method\n",
    "        * 'Wcycle'      W-cycle MultiGrid Method\n",
    "        * 'Fcycle'      Full cycle Multigrid Method\n",
    "        * 'cg'     mg-Preconditioned Conjugate Gradient\n",
    "        * 'minres' mg-Preconditioned Minimal Residual Method\n",
    "        * 'gmres'  mg-Preconditioned Generalized Minimal Residual Method\n",
    "        * 'bicg'   mg-Preconditioned BiConjugate Gradient Method\n",
    "        * 'bicgstable' mg-Preconditioned BiConjugate Gradient Stabilized Method\n",
    "        * 'bicgstable1' mg-Preconditioned BiConjugate Gradient Stabilized Method\n",
    "        The default setting is 'cg' which works well for SPD matrices. For\n",
    "        non-symmetric matrices, try 'gmres' and for symmetric but indefinite\n",
    "        matrices, try 'minres' or 'bicg' sequences.\n",
    "        The string option.solver is not case sensitive.\n",
    "    - option.preconditioner:  multilevel preconditioners including:\n",
    "        * 'V'   V-cycle MultiGrid used as a Preconditioner\n",
    "        * 'W'   W-cycle MultiGrid used as a Preconditioner\n",
    "        * 'F'   Full cycle Multigrid used as a Preconditioner\n",
    "        * 'bpx' BPX-Preconditioner\n",
    "    - option.printlevel: the level of screen print out\n",
    "        * 0: no output\n",
    "        * 1: name of solver and convergence information (step, err, time)\n",
    "        * 2: convergence history (err in each iteration step)\n",
    "\n",
    "\n",
    "`[x,info] = mg(A,b,elem)` also returns information of the solver\n",
    "\n",
    "    - info.flag:\n",
    "        * 0: mg converged to the desired tolerance tol within maxIt iterations\n",
    "        * 1: mg iterated maxIt times but did not converge.\n",
    "        * 2: direct solver\n",
    "    - info.itStep: the iteration number at which x was computed.\n",
    "    - info.time: the cpu time to get x\n",
    "    - info.err: the approximate relative error in the energy norm in err(:,1) and the relative residual norm(b-A*x)/norm(b) in err(:,2). If flag is 0, then max(err(end,:)) <= tol.\n",
    "    - info.stopErr: the error when iteration stops"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples: Jump Coefficients\n",
    "\n",
    "In [Introduction to Fast Solvers](solverintroduction.html), several examples has been presented for linear and quadratic elements. Here we present a harder example on a 3-D elliptic equation with jump coefficients. This documentation is based on `example/solver/Poisson3jumpmgrate.m`. Run this example to get more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ -\\nabla \\cdot (\\omega\\nabla u) = f\\quad  \\text{ in } \\Omega=(-1,1)^3$$ \n",
    "$$u = 1 \\text{ on } x=1, \\qquad u=0  \\text{ on } x=-1$$ \n",
    "$$\\omega\\nabla u \\cdot n = 0$$ on other boundary faces.\n",
    "\n",
    "The diffusion coefficent $\\omega$ is piecewise constant with large jump:\n",
    "- $\\omega(x) = 1$ if $x\\in (-0.5, 0)^3$ or $x\\in (0,0.5)^3$ and \n",
    "- $\\omega = \\epsilon$ otherwise.  \n",
    "\n",
    "![Domain](jump3d.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Reference**\n",
    "\n",
    "Xu, J, and Y Zhu. “Uniform Convergent Multigrid Methods for Elliptic Problems with Strongly Discontinuous Coefficients.” M3AS 18, no. 1 (2008): 77–106."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Setting\n",
    "% mesh\n",
    "[node,elem,HB] = cubemesh([-1,1,-1,1,-1,1],1);\n",
    "bdFlag = setboundary3(node,elem,'Dirichlet','(x==1) | (x==-1)');\n",
    "mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);\n",
    "% option\n",
    "option.L0 = 3;\n",
    "option.maxIt = 4;\n",
    "option.elemType = 'P1';\n",
    "option.printlevel = 1;\n",
    "option.plotflag = 0;\n",
    "option.dispflag = 0;\n",
    "option.rateflag = 0;\n",
    "% pde\n",
    "pde = jumpmgdata2;\n",
    "global epsilon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 12,   err = 8.24e-09,   time = 0.13 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 13,   err = 2.30e-09,   time = 0.21 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 13,   err = 4.14e-09,   time =    2 s\n",
      "\n",
      " Table: Solver MGCG for epislon = 1.00e-01 \n",
      " #Dof   Steps Time      Error    \n",
      "\n",
      "  4913   12   0.13   8.2382e-09\n",
      " 35937   13   0.21   2.2998e-09\n",
      "274625   13      2   4.1353e-09\n",
      "\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 13,   err = 9.52e-09,   time = 0.06 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 14,   err = 7.92e-09,   time = 0.25 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 15,   err = 4.04e-09,   time =  1.8 s\n",
      "\n",
      " Table: Solver MGCG for epislon = 1.00e-02 \n",
      " #Dof   Steps Time      Error    \n",
      "\n",
      "  4913   13   0.06   9.5205e-09\n",
      " 35937   14   0.25   7.9168e-09\n",
      "274625   15    1.8   4.0409e-09\n",
      "\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 14,   err = 2.22e-09,   time = 0.03 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 15,   err = 7.20e-09,   time = 0.22 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 16,   err = 3.47e-09,   time =  1.9 s\n",
      "\n",
      " Table: Solver MGCG for epislon = 1.00e-03 \n",
      " #Dof   Steps Time      Error    \n",
      "\n",
      "  4913   14   0.03   2.2209e-09\n",
      " 35937   15   0.22   7.1950e-09\n",
      "274625   16    1.9   3.4668e-09\n",
      "\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 14,   err = 2.37e-09,   time = 0.04 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 16,   err = 4.41e-09,   time = 0.22 s\n",
      "Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 17,   err = 4.08e-09,   time =    2 s\n",
      "\n",
      " Table: Solver MGCG for epislon = 1.00e-04 \n",
      " #Dof   Steps Time      Error    \n",
      "\n",
      "  4913   14   0.04   2.3686e-09\n",
      " 35937   16   0.22   4.4141e-09\n",
      "274625   17      2   4.0804e-09\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%% MGCG (multigrid preconditioned CG) solver\n",
    "for k = 1:4\n",
    "    epsilon = 10^(-k);\n",
    "    [err,time,solver,eqn] = femPoisson3(mesh,pde,option);\n",
    "    fprintf('\\n Table: Solver MGCG for epislon = %0.2e \\n',epsilon);\n",
    "    colname = {'#Dof','Steps','Time','Error'};\n",
    "    disptable(colname,solver.N,[],solver.itStep,'%2.0u',solver.time,'%4.2g',...\n",
    "                      solver.stopErr,'%0.4e');\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multigrid Vcycle Iteration \n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 31,   err = 7.92e-09,   time = 0.14 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 34,   err = 8.42e-09,   time = 0.33 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 35,   err = 8.12e-09,   time =  2.8 s\n",
      "\n",
      " Table: Solver V-cycle for epislon = 1.00e-01 \n",
      " #Dof   Steps Time      Error    \n",
      "\n",
      "  4913   31   0.14   7.9248e-09\n",
      " 35937   34   0.33   8.4179e-09\n",
      "274625   35    2.8   8.1222e-09\n",
      "\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 83,   err = 8.29e-09,   time = 0.09 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 122,   err = 8.88e-09,   time = 0.92 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 151,   err = 9.18e-09,   time =  9.6 s\n",
      "\n",
      " Table: Solver V-cycle for epislon = 1.00e-02 \n",
      " #Dof   Steps  Time      Error    \n",
      "\n",
      "  4913    83   0.09   8.2928e-09\n",
      " 35937   122   0.92   8.8789e-09\n",
      "274625   151    9.6   9.1830e-09\n",
      "\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 122,   err = 9.89e-09,   time = 0.23 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 200,   err = 2.44e-07,   time =  1.3 s\n",
      "NOTE: the iterative method does not converge! \n",
      "Multigrid Vcycle Iteration \n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 200,   err = 4.59e-05,   time =   12 s\n",
      "NOTE: the iterative method does not converge! \n",
      "\n",
      " Table: Solver V-cycle for epislon = 1.00e-03 \n",
      " #Dof   Steps  Time      Error    \n",
      "\n",
      "  4913   122   0.23   9.8922e-09\n",
      " 35937   200    1.3   2.4413e-07\n",
      "274625   200     12   4.5888e-05\n",
      "\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:     4913,  #nnz:    28747, smoothing: (1,1), iter: 130,   err = 8.91e-09,   time = 0.37 s\n",
      "Multigrid Vcycle Iteration \n",
      "#dof:    35937,  #nnz:   230043, smoothing: (1,1), iter: 200,   err = 1.16e-06,   time =  1.3 s\n",
      "NOTE: the iterative method does not converge! \n",
      "Multigrid Vcycle Iteration \n",
      "#dof:   274625,  #nnz:  1838395, smoothing: (1,1), iter: 200,   err = 2.17e-04,   time =   11 s\n",
      "NOTE: the iterative method does not converge! \n",
      "\n",
      " Table: Solver V-cycle for epislon = 1.00e-04 \n",
      " #Dof   Steps  Time      Error    \n",
      "\n",
      "  4913   130   0.37   8.9149e-09\n",
      " 35937   200    1.3   1.1554e-06\n",
      "274625   200     11   2.1676e-04\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%% V-cycle solver\n",
    "for k = 1:4\n",
    "    epsilon = 10^(-k);\n",
    "    option.mgoption.solver = 'Vcycle';\n",
    "    [err,time,solver,eqn] = femPoisson3(mesh,pde,option);\n",
    "    fprintf('\\n Table: Solver V-cycle for epislon = %0.2e \\n',epsilon);\n",
    "    colname = {'#Dof','Steps','Time','Error'};\n",
    "    disptable(colname,solver.N,[],solver.itStep,'%2.0u',solver.time,'%4.2g',...\n",
    "                      solver.stopErr,'%0.4e');    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "For 3-D jump coefficients problem, MGCG (multigrid preconditioned CG) solver converges uniformly both to the mesh size and the ratio of the jump. \n",
    "\n",
    "V-cycle alone doesn't converge for small epsilon and small h."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "120px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
